Optimal Filling of Shapes 
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We present filling as a type of spatial subdivision problem similar to covering and packing. Filling 
addresses the optimal placement of overlapping objects lying entirely inside an arbitrary shape so 
as to cover the most interior volume. In n-dimensional space, if the objects are polydisperse n-balls, 
we show that solutions correspond to sets of maximal n-balls. For polygons, we provide a heuristic 
for finding solutions of maximal discs. We consider the properties of ideal distributions of A'^ discs 
as N ^ oo. We note an analogy with energy landscapes. 



Packings of non-overlapping objects such as monodis- 
perse or polydisperse spheres, ellipsoids, or polyhedra 
have been long-studied by physicists and mathematicians 
[1-11]. Coverings of shapes by overlapping objects are 
also of interest in many physical settings[12-16]. Pack- 
ing and covering are two familiar examples of problems in 
the subdivision of space subject to prescribed constraints. 
Whereas in the packing problem, objects packed in a 
given shape are not allowed to overlap each other or the 
shape boundary, in the covering problem objects overlap 
both each other and the shape boundary in an effort to 
maximally cover a given shape. In both problems, the 
objects may be monodisperse or polydisperse in size. In 
covering, the objects are typically n-balls (discs, in 2D); 
in packing, the objects may be of any shape. 

In this Letter, we present a new type of spatial subdi- 
vision problem which can be viewed as intermediate be- 
tween the packing and covering problems. We define fill- 
ing as the problem of packing overlapping objects inside 
of a defined shape so as to cover the interior volume with- 
out extending beyond the boundary of the shape (Fig. 1). 
We are primarily interested in the optimal filling of an n- 
dimensional shape, characterized by a well-defined n — 1 
dimensional surface, with N polydisperse n-balls. Specif- 
ically, we seek the optimal placement and radii of the 
n-balls for a given N. 

The notion of filling arises from the problem of model- 
ing anisotropic nanoparticles as rigid bodies composed 
of a sum of isotropic volume-excluding potentials [17]. 
Other applications are the problem of irradiating a tumor 
with the fewest number of beam shots, while controlling 
the beam diameter, but without damaging surrounding 
tissue [18]; using time-delayed sources to create shaped 
wavefronts; combining precision-placed explosives with 
tunable blast radii; positioning proximity sensors with 
defined radii; cell phone and wireless network coverage; 
or any problem of ablation or deposition where one has a 
sharp impenetrable boundary and a radially tunable tool. 
It also can be related to the coarsening (due to Ostwald 



ripening) of wet foams packed in containers with non- 
wetting surfaces. 
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FIG. 1. The problem of filling a shape, such as a triangle, 
may be viewed as intermediate between the familiar packing 
and covering problems. 

Here we show that, to optimally fill an n-dimensional 
shape with N n-balls of varying radii, only solutions of 
maximal n-balls (that is, balls whose centers lie on the 
medial axis of the shape) need be considered. It will fol- 
low that the dimension of the solution space is reduced 
from n + 1 to n — 1 . Secondly, we consider optimal fillings 
of polygons with N polydisperse overlapping discs. For 
this two-dimensonal case, we present a heuristic for the 
numerical generation of optimal solutions for arbitrary 
N. Thirdly, by considering how discs are optimally dis- 
tributed in polygons as iV — cxi, we derive exact analyt- 
ical expressions for the spatial distribution of the discs 
and show that the fraction of unfilled area for optimal 
solutions vanishes like 1/7V^. The analytical expressions 
may be used to approximate solutions for finite but large 
N. We derive an exact expression for the fractional al- 
location of discs over the three medial axis branches of 
a triangle as iV oo. We discuss how solutions for 
n — 2 provide insight into the filling problem of gener- 
alized shapes in arbitrary dimensions. We also note an 
interesting connection to energy landscapes. 
Reducing the dimension of the solution space using max- 
imal n-balls. For optimal filling solutions, the objective 
function to be maximized is defined to be the volume of 
the union of a set of N n-balls constrained to the interior 
of a shape G. The upper bound of this function is the to- 
tal volume of G. The optimal A'' = 1 filling is the largest 
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n-ball that can fit in G. Given a set of TV n-balls, the 
contribution of a single bah to the total filling is equal to 
the volume of the ball minus the fractional share of the 
volume of any overlap with the A^ — 1 other balls (SI, §1). 
To find optimal solutions for G, we find it is not necessary 
to consider the space of all possible n-balls contained in 
G. Importantly, we show that only the maximal n-balls 
need be considered. 

First introduced by Blum[19, 20] as a "topologi- 
cal skeleton", the medial axis is a reduction of an 
n-dimensional shape into an n — 1-dimensional space, 
M(G), the locus of centers of the maximal n-balls. A 
maximal n-ball is an n-ball completely contained in the 
shape tangent to the shape boundary at two or more 
points. Also, a maximal n-ball is a ball contained com- 
pletely in G but not contained in any other ball in G. 

The radius function r is a continuous, non-negative 
function defined at each point of M{G) as the radius of 
the maximal n-ball centered at that point. The medial 
axis and the radius function comprise a complete shape 
descriptor[20] and can be used to reconstruct G. Every 
maximal n-ball in G can be represented as a unique point, 
its center, on M{G). If each n-ball in a proposed opti- 
mal filling is replaced by a maximal n-ball containing it, 
the new solution will fill an equivalent or greater area of 
G. Therefore, only solutions of maximal n-balls need be 
generated. Finding solutions is reduced to finding center 
points on an n — 1-dimensional hypersurface. 
Heuristic for 2D shapes. For a planar shape G, M{G) 
is the one-dimensional planar graph that is the locus of 
the centers of the maximal discs (2-balls) of the shape. 
M(G) is a set of 1-manifolds, or branches, plus connect- 
ing branch points and terminating end points [20]. For a 
convex polygon, M(G) is composed only of straight seg- 
ments. For polygons[21], M{G) is composed of straight 
segments and, possibly, parabolic curves. Various algo- 
rithms exist to compute the medial axes of convex or 
concave polygons [22, 23]. The medial axes of a pentagon 
and a concave polygon are shown in Fig. 2. The filled 
area of all the N — \ maximal disc fillings is shown from 
the side in Fig. 2 for both shapes. The global maximum 
is synonymous with the largest disc inscribable in G. In 
a convex shape, there is only one maximum. In a concave 
shape, there can be many. 

The neighbors of a disc A in M{G) are the set of centers 
that can he reached along any path in M{G) originating 
at the center of A without traversing another center. In 
2D, we find that the change in the total filling due to 
locally displacing a maximal disc center on M(G) is a 
function of only the change in the overlaps of the maximal 
disc with its neighbors(SI, §1). In Fig. 2, the centers of 
the neighbors of the largest disc of the concave polygon 
are enclosed by dashed boxes. 

Traps are a special set of points on the medial axis, 
which are important in optimal filling solutions because 
they are often occupied by centers. To see this, we con- 
sider the behavior of the objective function around a trap. 
We begin by observing that the filling function is piece- 



wise first-order continuous, with fixed points of first-order 
discontinuity that we refer to as junctions. These are 
points where the radius function and/or path in M{G) is 
discontinuous. In 2D all branch points are junctions. If 
a junction is a local maximum with respect to moving a 
single disc center (all other disc centers held fixed), then 
the junction is a trap, because small displacements in 
the centers of discs neighboring the discontinuity do not 
displace the position of the local maximum. In Fig. 2, 
the centers caught in traps in the optimal filling solu- 
tions are shown as red open dots. We observe that disc 
centers in traps tend to be common features in optimal 
filling solutions, and even a fixed feature when N is large. 

We now propose a solution strategy whereby disc cen- 
ters are distributed onto M{G) and the local filling max- 
imum is found by a gradient method. If the maxima of 
enough distribution samplings are generated such that 
all the local maxima can be enumerated, then the global 
filling maximum is among them. We propose the fol- 
lowing heuristic that accomplishes this efficiently while 
also greedily using the N — 1 solution to intelligently 
reduce the number of distributions to be searched: (1) 
The medial axis is divided into K pieces, branches with 
monotonically increasing radius functions and the junc- 
tions connecting them. (2) It is assumed that there is 
at most one local maximum per way, W, of partitioning 
the N discs over the K pieces (branches and junctions), 
W = {ni}f^ where N = J^f f^i and Ui is the number of 
discs on the ith piece, n^ G N. (3) It is further assumed 
that, given the optimal way of partitioning N — 1 discs, 
{n-}f- , the optimal way of partitioning A'' discs is nearby, 
where nearby means | n^ — n^ | is small, and that if 
the discs assigned to a given piece are decreased in num- 
ber, the pieces that have discs increased in number have 
a minimal distance (counted by number of connecting 
pieces) to tlic decreased piece. (4) The local maxima of 
the nearby ways are generated using any local maximum 
finding technique (e.g. active set or sequential quadratic 
programming optimization schemes [24]). The best local 
maximum found is presumed to be the optimal N filling 
solution for the shape. Note that during the local maxi- 
mum search, the searched space is first-order continuous 
because the fixed points of first-order discontinuities are 
in the set of junctions. Non-fixed points of first-order 
discontinuity that correspond to the point of tangency of 
discs in G are not generally relevant because they cannot 
be points of local maximum. 

This heuristic is made more efficient by taking advan- 
tage of center-occupied traps and the dependence of the 
filling function on the nearest neighbors. When a trap 
is occupied by a center in the TV — 1 solution, the phase 
space of centers can be divided into independent sub- 
spaces. If it is known (or guessed) that the best solution 
for A'^ also includes a center in the trap, then the sub- 
parts of M(G) connected only by the center-occupied 
trap can be searched independently. Rearrangements of 
centers in one sub-part cannot affect the best arrange- 
ment of centers in another if they are connected only by 




FIG. 2. Two polygons and their medial axes (green line segments): (a) pentagon and (b) concave polygon. Below each is a side 
view of the landscape of the N = 1 total filled area function (the area of both shapes is normalized to 1). For the pentagon, 
the horizontal extents of the shape are shown by dashed lines. For the concave shape, the local filling maxima of the A'^ = 1 
solution are shown by dashed lines. To the right of each shape is an optimal filling with 21 discs. Disc centers in traps are 
shown as red open dots. For the concave polygon, dashed squares enclose the centers of neighbors of the largest maximal disc. 



the center-occupied trap. 

We implemented this heuristic for polygons, which 
have only a few types of medial axis branches. The sets of 
discs shown in Figs. 2 and 3 are the best solutions found 
by the heuristic and verified by a genetic algorithm [2 5]. 
The heuristic generally produces a superior solution to 
the genetic algorithm, as long as a large enough neigh- 
borhood of ways is considered. One surprising finding is 
that the largest disc that fits in G, i.e. the = 1 solu- 
tion, is not always part of the solution for iV > 1. In Fig. 
3, for example, the N ^ 5 solution for the trapezoid does 
not include the A^=:l solution. 

Optimally filling a polygon as N oo. It is instructive to 
examine how the optimal filling of a shape converges to 
the total volume of the shape as — oo. We find that 
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FIG. 3. Maximal filling solutions for polygons 



this limit can be solved exactly for polygons. The M{G) 
of a polygon can be divided into branches of only three 
types, (1) straight segments with linear radius functions, 

(2) parabolic curves with quadratic radius functions, and 

(3) straight segments with square root radius functions 
(SI, §2). Immediately, case (3) can be ignored because 
no disc center on such a curve fills more area than that 
filled by two discs placed at the ends of the curve. In 
ideal solutions. Case (3) type curves are empty except 
for the ends of the curve. Thus we need only derive an 
expression for the relative density of centers as a function 
of position on the other two branch types. Along each 
branch type, the distribution of optimal filling solutions 
approaches this expression as iV — >■ oo. 

Let p{t) represent the density of centers, r{t) the radius 
function, K{t) the local curvature, and {x{t),y(t)) the 
position along a parameterized branch of M(G'), with 
i G [ta, tb] ■ Given an expression for the unfilled area Ai 
along the path i of the form 



A, 



(1) 



where C'i is a function to be determined, we wish to de- 
termine the function p(t) that minimizes this area con- 
strained by = pdt. Note that if we sum the un- 
filled areas Ai over all of M{G), then the filled area is 
Ag — where Ac is the area of G. This variational 

problem can be solved by constructing the Lagrangian 



Qit)j.~\p 



dt 



(2) 



and taking the pointwise derivative with respect to p{t), 



dp{r) 



It! + lA-.C^it) - a) S{t - r)dt = 0. 



(3) 



This equation is solved by functions p that satisfy 

-2ait)+pi.ait)-p^\ = (4) 
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Note that p = [j^^y^^ satisfies this equation. 

For Case (1), where y{t)) = At + B, and r{t) = 

ct + ro, for constants c, ro > 0, where ta > 0, we derive 
(SI, §3) 



3/2 



(5) 



It follows that yO cx r . 

For Case (2), where {x{t),y{t)) = {2rot,rot'^), r = 
ro(l+f^), where ro is the minimum of the radius function, 
and K{t) = 2^ (1 + ^^^.^^ ^gj^ 



It follows that p oc (l + t^ 
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(6) 
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For both Case (1) and Case (2), the distribution of 
centers follows a power law with respect to the local ra- 
dius function. Wc observe that centers on M (G) will 
be distributed more densely where the radius function is 
smaller. Given p = por~", for a = 1/3 or 5/6, po can be 
determined from 



PO =N (//; r-^dt) ' = N/Ro. 



(7) 



Ro is then a constant determined by the radius function 
of the branch. 

For Case (1), we observe that the distribution of cen- 
ters on the medial axis branch is also scale-free. The 
distribution of centers follows a power law with respect 
to the distance from the vertex (where t = 0) of the 
polygon. Eq. 1 becomes 
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J'" RlC{K,r',r)dt 



(8) 



where C is the evaluated integral. Thus, as ^ oc for a 
system of ideally distributed centers, the filling converges 
to the area of G with an asymptotic error proportional 
to iV~^. It can be presumed that all shapes that can be 
approximated by polygons with an increasing number of 
sides also converge with an N~'^ error term. 

If we divide M{G) into k branches we can predict the 
partitioning of the discs over the branches as A/' — >■ oo. 
The fraction of discs on a given branch i is (see SI, §5) 



{Ci 



\l/3 



(Cl)l/3 + (C2)V3 + ... + (C,)l/3- 



(9) 



For a triangle, which is composed of three Case (1) 
branches, the fraction of the discs on a given branch can 
be solved analytically: 



cot (6*1) + cot (^2) + cot (6*3) 



(10) 



where 61 , 62 and 63 are the internal angles of the triangle, 
each of which is associated with a branch of the medial 
axis. From Eq. 10, it is clear that the optimal solution 
preferentially populates medial axis branches associated 
with smaller internal angles. This can be observed in the 
triangles in column one of Fig. 3. 

/discussion. A convenient framework for visualizing fill- 
ing solutions of a 2Z) shape is to consider the unfilled area 
as the "energy" of the system, and the force acting on a 
disc center as the negative gradient of this energy. The 
force acting on a single center can be divided into two 
parts, a force due only to the local radius function, and a 
purely repulsive short range force (also dependent on the 
radius function) between a center and its neighbors (SI, 
§1C). The range is defined by where two discs overlap. 
Because this energy function is not first-order continu- 
ous, these force definitions have discontinuities. A center 
in a trap, therefore, has a restoring force in each path di- 
rection away from the trap. A center not in a trap is at a 
point where the local forces balance. The filled area plots 
of Fig. 2, therefore, is like an inverted energy landscape 
of the A^ = 1 system of the two polygons. 

In this work we defined the filling problem for arbi- 
trary dimensions. Although we examined the details of 
its solution structure only in two spatial dimensions, we 
can extrapolate some of our findings to higher dimen- 
sions. For instance, in the case of polygons, we found 
that first order continuous manifolds meet at lower di- 
mension manifolds where centers are trapped. It is natu- 
ral to expect similar behavior for D > 2. Further, higher 
dimensional polyhedron shapes are also likely to have 
manifolds with scale- free solutions. Such behavior will 
be investigated in future publications. 
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